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Abstract 

We explore a new framework for describing the kinetics of a heterogeneous chemical 
reaction where two particles of the same chemical species form a reaction product of another 
chemical species on the surface of a seed particle. Traditional treatments neglect the effect 
of statistical fluctuations in populations. We employ techniques in a manner analogous to 
the treatment of quantum systems to develop a stochastic description of processes beyond 
the mean field approximation. 

1 Introduction 

Traditional models employing the evolution of the mean population of species in a system provide 
a good enough description of various physical processes such as heterogeneous chemical reactions 
and heterogeneous nucleation of aerosols as long as the average particle number is large. However, 
when the mean populations involved are small we expect to see a significant deviation from the 
solution to the classical equations. 

As an example consider the chemical reaction of two hydrogen atoms on the surface of a dust 
particle in interstellar space [13j . Atoms can be adsorbed onto and evaporated from the grain 
particle. The adsorbed molecules will diffuse on the surface of the seed and eventually collide 
with another reactant to form diatomic hydrogen. Under interstellar conditions the rate of the 
adsorption of atoms will be small compared to the reaction rate. If the mean population of hy- 
drogen atoms is small, statistical fluctuations brought about by gain and loss, and by the random 
diffusion of atoms on the surface of the dust particle, are important. Yet the traditional model 
does not include the treatment of such fluctuations. Several attempts have been introduced to 
resolve that problem; among them the Monte Carlo approach, the modified rate approach, the 
direct master equation approach and the Gauge Poisson representation approach — see [21 and 
. Related studies have been carried out in [H [11] . Analytical solutions to a master equation 
approach have been found for the steady state case[Tl [HI [H]. Our work is a first step to an al- 
ternative approach for describing the kinetics of a heterogeneous chemical reaction which can be 
extended to other areas where a similar problem occurs, for example when computing the rate of 
nucleation processes taking place on small particles. Again, the number of adsorbed molecules is 
small and statistical fluctuations need to be taken into account. In [16j the steady state solution 
was studied employing a description using master equations instead of mean population dynam- 
ics. A nucleation rate lower than the classical equations predicted was obtained. Therefore, it is 
important to develop a model by which we replace the classical equations with a set of stochastic 
equations. 
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For the correct treatment of population fluctuations we employ methods of Quantum Field The- 
ory [SI [S] . Starting from a master equation we introduce a spatial lattice where the microstates 
of the system correspond to a set of occupation numbers at each lattice site. A Fock space 
is constructed using annihilation and creation operators at each lattice site. By means of this 
set-up it is easy to show that the master equation is equivalent to a Schrodinger equation with 
imaginary time. This enables us to employ techniques originally developed in order to describe a 
quantum mechanical system where the fluctuations are due to a quantum uncertainty. We obtain 
the average particle population of the classical many-body system by developing a mechanism for 
computing expectation values of observables analogous to Feynman's path integral formulation 
[21 [18] . In this paper we will not concern ourselves with renormalisation group analysis although 
studies in that direction have been undertaken by various people, for example [121 [13 [21] ■ Intro- 
ducing a stochastic variable [20], a mathematical trick — for a nice review paper we refer to [21] — 
helps with the evaluation of the expression for the expectation values. The complex fluctuating 
solutions to a set of constraint equations, which are stochastic partial differential equations, are 
then averaged over all realisations of the stochastic noise. For numerical investigations, the solu- 
tions to the constraint equations can be generated by various numerical schemes |14j . The path 
integral average is computed using Monte Carlo methods [17]. The Code is written in C and com- 
putations do not take longer than a few seconds up to a few minutes running on a standard laptop. 

In the sequel we will concentrate on a heterogeneous chemical reaction where two reaction part- 
ners of the same sort of particle type react on the surface of a seed particle to form a reaction 
product. 

For readers more interested in the physical content than in the mathematical details of the 
formalism we recommend to skip the first few pages and instead have a look at the most im- 
portant equations on page 8 and continue from there. Equation (|42|) which is a path integral 
average (PIA) gives the average particle density of the various chemical species once the complex, 
fluctuating solutions to the constraint equations ([40]) and (j41|) have been found and inserted into 
the PIA. Together with the correlations for the stochastic noise occuring in the PIA this 
forms a complete set of equations. 



2 Mathematical Techniques 

2.1 Master Equation And Schrodinger-like Equation 

We concentrate on chemical reactions of type A + A — > C, that is situations in which the 
atoms adsorb onto grain particles where they can associate with themselves to produce diatomic 
molecules. If the number of reactive species on an individual grain is small, traditional rate 
equations will fail to accurately describe the diffusive chemistry occuring on the surface of the 
grain particle. We start our investigations by determining a master equation that describes such 
a heterogeneous chemical process. 



A general master equation can be written as 

= J2 Tn^rnP{n) - ^ r„^„P(m), (1) 
n n 

where Tn^m represents the transition amplitude or propagator from a microstate n to a microstate 
m and P{m) is the probability to find the system in state m. Considering a d-dimensional lattice L 
with lattice constant I, the microstates correspond to the occupation numbers {Ni} — {Ni, N2, ...} 
at each lattice site i. 



The chemical reaction of pairs of species A to form a product of species C on a particle or 
droplet , A + A — > C, is modelled by the following master equation 
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+ ^ (P({iVA}..., A^c. - 1, i) - {^c}; t)) 

i 

- iVA.(A^A. - l)-P({A^A},{iVc};i)) 
+ XaY. {{Na, + 1)P{...,Na, + 1, {Nch t) - Na,P{{Na}, {NcY t)) 

i 

+ XcYl ((^c. + + 1, ■■■■,t) - NcPHNa}, {Nc};t)) 

i 

+ DAj2{i^A, + ^)P(-'^A, + ^^^A,-l,-,{Ncht)~NA,P{{NA}{Nc};t) 

(u> 

+ {Na, + 1)P{...,Na, - 1,Na, + 1, {Nc};t) - Na,P{{Na}, {Nc};t)) 

+ DcY, ((^c. + 1)P{{Na}, Nc, + 1, Nc, - 1, t) - Nc,P[{Na}{NcY t) 

+ [Nc, + 1)P{{Na}, -.Nc, - 1, Nc, + 1, t) - Nc,P{{Na}, {Nc}; t)), 

(2) 

This equation describes the evolution of the probabihty distribution P{{NA},{Nc};t) for the 
total number of adsorbed molecules {Na} of species A and for the number of reaction products 
{Nc} of species C. The symbols NAi,Ci denote the numbers of A or C molecules at lattice site i, 
respectively. The constants jA,c, ^ and Xa,c are rate coefficients, Da,c is the diffusion constant 
and V stands for the volume of the droplet. The rate coefficient jA,c is called source rate, the 
rate coefficient Xa,c is called evaporation rate and n is known as the reaction rate. 

In our model, the chemical reaction is taking place on a d-dimensional lattice, allowing for mul- 
tiple occupancy on each site. This configuration is also called bosonic representation. 
The changes in population which we consider are caused by: 

(a) absorption of a molecule of species A from outside the grain particle (first line in the above 
equation), and absorption of a molecule of species C from outside the grain (second line in 
the above equation), 

(b) binary reaction on the surface of the grain (third and fourth line in the above equation), 

(c) evaporation of a molecule of species A from the grain (fifth line in the above equation), 
and evaporation of a molecule of species C from the grain particle (sixth line in the above 
equation) , 

(d) particle hopping of a molecule of species A from site i to site j (seventh line in the above 
equation) , and particle hopping of a molecule of species C from site i to site j (eighth line 
in the above equation), 

(e) particle hopping of a molecule of species A from site j to site i (ninth line in the above 
equation), and particle hopping of a molecule of species C from site j to site i (last line in 
the above equation). 

The summation in the fourth and fifth line of the above equation ^ is taken over nearest 
neighbour sites only. The factors {Na,c + 2), {Na,c + ^),Na,c, {Na,c — 1) describe the number 
of ways of choosing the particles involved in the considered process. In the continuum limit the 
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particle hopping from one site to another corresponds to the diffusion of the particles. 
The initial condition is chosen corresponding to a Poissonian distribution on each site 

P{{NAUNcht^ 0) = e--(°)--(°) n "^^"Cl^c^r^' ' 

i ^ ^ 

where ^^^^(O) is the average occupation number per lattice site for the A or C particles respec- 
tively. 

In the next step we will apply the methods of second quantisation [51 [5] . We will rewrite the 
master equation as a Schrodingcr-likc equation for a many-body wave function. This approach 
can be justified by noting that, first of all, the master equation is a differential equation of first 
order with respect to time. The second reason to suggest the treatment of the master equation 
according to the second quantisation is that the master equation is linear in the probability. 

In order to simplify the notation we will suppress the dependence on space coordinates x = 

{xi,X2, ...,a;d). We will be working in an appropriate space, the Fock space. A Fock space 

is a Hilbert space made from the direct sum of tensor products of single-particle Hilbert spaces 

n 

OC 

^.(^) =0^.^®", (4) 

n=Q 

with Si, a symmetrising (for the case of bosons) or antisymmetrising (for the case of fermions) 
operator. The Fock space is constructed by introducing the following operators at each lattice 
site i 

a^, i G L : creation operator, 
aj, i G L : annihilation operator. 



which satisfy the commutation relationships 

^[aj, aj] ^(a.a^- - a^a^) = Sij. (5) 
The vacuum state |{0}) is defined by 

a. |{0}) = |{0}) V*eL, (6) 

with 

|{0}):=(g)|0,), (7) 

j 

where |0)j denotes the vacuum state in a single-particle Hilbert space. 

The master equation ([2]) is equivalent to the Schrodinger-like equation — a Schrodinger equa- 
tion with imaginary time — 

^|*)a+a-c = -HA+A^c[aA,, a^j, acfc, acJ|*)A+A-.c (8) 
with the many-body wave function 

\^)a+a^c ■■= E Pi{NA},{Ncht)ll{tA.f^'itc.f^^m), (9) 

{Na},{Nc} i 
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and the Hamiltonian operator 



Me{A,C} i 

,2 2 



K 

V 



II ( ac. - a^^ ) a. 



+ X! -Dji/ ^(aA/, - aMj)(aM. - ajv/J. (10) 

Me{A,C} (ij) 

For verification of tlie above statement one lias to insert the states Yiii^Ai)^^'- (^Ci)^'^'^\{0}) 
on both sides of the master equation — equation ^ — and sum over the set of all occupation 
numbers {Na} and {Nc}- In general, the time-evolution operator H is not necessarily Hermitian. 
For the timebeing let us suppress the subindices that identify the particle type. The form of the 
many-body wave function ^ can be made plausible when considering the state vector \Ni) at 
site i e L, namely 

1^^) =ar |0). (11) 

It holds that 

a, lA^,) = A^.liV, -1), a, l^'*) = l^. + l>- (12) 

2.2 Expectation Values Of Observables 

We are interested in obtaining the expectation values for various observables, especially the 
average number density. The expectation values of observables D are given by 

(0):=5]0({iV,;})P({7Vj;i). (13) 

We want the expectation values of the observables — diagonal in the occupation number basis — 
to be linear in the probabilities. After some straightforward computation — see for example [23j — 
one can see that the above expression is equivalent to 

(D) = ({P}|D|M^(i)), (14) 

where 

({P}| := ({0}|e^3^^ : projection state. (15) 
The projection state obeys the relation 

({P}|{0}) = 1. (16) 
By definition, the projection state is a left eigenstate of all creation operators with unit eigenvalue 

({P}|t,= ({P}|, VzeL. (17) 
Furthermore, ({P}|'I'(t)) = 1. Conservation of probability requires that ({P}|H — ({0}|. 

We break the time interval [tajtr] into T short slices of duration At = . At each time 

slice we insert a complete set of coherent states — see, for example, pj. Coherent states, \Ci{t)), 
are right eigenstates of the annihilation operator 

a, |C,(i)) = $,(t)|C,(i)), ieh (18) 

where the eigenvalue is a complex function. The duals (C; (t) \ are left eigenstates of the creation 
operator 

(C,(i)| a,^ (C,(<)|$*(<), ieL. (19) 
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We have 



(C,(t)| := (0|e-^l*»*(*)l'+**(*)"'. (20) 
The coherent states are over-complete. Stih, we can use them to create the identity 

1 = i / d[Rei<P,)]d[Im{<P,)]\C,it)){C,{t)\, (21) 



for a single lattice site i e L, and for multiple lattice sites accordingly 

1 = / n {^d[Remd[Imm) l{C})({C}|, (22) 

with |{C}) = |Cj). Let us recall the formula for the expectation values of observables 

(0) = ({P}|0|vl/(i)) = ({0}|e^-^-Oe-''*|vI/(0)), (23) 
where the initial many-body wave function takes the form — see equations ([3|) and ([9]) — 

l^'(O)) := e"("K2:,i»-i)|{0}) (24) 
We observe the following proportionalities 

({P}|(x({C(t)}k..i = ({0}|e-HE.a.^ 

\^{t^Q)) (X |{C(t)})$^=„(o) =e-^l"(°)l'+"(")^-*M{0}), (25) 

for all admissible values of j. Therefore, one can recast the equation for the expectation values 
([23)1 into 

(0) « ({Ci(t)}|Oe-«*|{C„(o)(i)}), (26) 

where 

({Ci(t)}| := ({C(t)}U..i, 

|{C„(o)(^)}) := |{C(t)})$,=„(o), (27) 

for all admissible values of j. According to the breakage of the time interval into time slices of 
small duration we rewrite the expression 

g-Ht ^ g-HAtg-HAt (T times) (28) 

occuring in the equation for the expectation values (pS)) and insert the identity as defined in ([^ 
between each factor. Then the discrete version of the expectation values of operators reads 

{0)d^screte'X J (j]_d[Re{<^^,r)]d[Im{<S>,,^)]j {{Ci}\0\{Cr}) X . . . 

■■■■ X ( n ({C Jle-^'^^KC-A.})) ({C.=o}|{C„(o)}>, (29) 

where we have labeled each time slice by a time index r S [0, At, 2 At, ■■■,T]. The normalisation 
constant has to be determined lateron. The consideration of the lattice expectation value is not 
sufficient if one is interested in long wavelength properties. In the formal continuum limit, we 
obtain 

\0) continous = (0) = hni (0) discrete • (30) 
At ►O 
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Next, we expand the exponential function for small At, neglect higher order terms in At and re- 
cast the continuous average. When interested in inclusive probabilities — e.g. the average number 
of particles at a given lattice site irrespective of the number of particles elsewhere — it is conve- 
nient to commute the factor of e^^^' through the operators D and H in (0). This has the effect of 

shifting a^a +1 using a= (a +l)e^. The operators are now normal ordered. It is valid that 
the operator and its normal ordered counterpart have the same expectation value if all creation 
operators occuring in the normal ordered operator are replaced by unity — see for example [23j . 

In particular, the density operator aa reduces to the annihilation operator a. Therefore, in the 
continuum limit the average particle density of the A particles is given by the stochastic average 
of the complex eigenvalues of the coherent state vectors under the annihilation operator, that is 
one chooses the operator to be ^a{x, t). 

We take the continuum limit — the dimensions of the constants are chosen by looking at the dis- 
crete Hamiltonian operator — via J2i — ^ / d^^i ^A,,Cj (t) — > *A,c(x, t)l° , ^*Ai,Cj (*) ~^ 

,c, 3A,c JAfi'l'^ and finally n^,c(0) — > 
nA.c((^)l'^ , where the newly introduced constants have the following dimension properties [^a.c{^i t)] 
m^^, [Da,c] = m^s"^ , [k] = m^s~^, [jA,c] = m~^s~^ and [71^^(7(0)] = in Standard Inter- 

national units. The objects ^\Q{x,t) and Xa,c are dimensionless. Notice that now ^A,c{^,t) 
scales like a density. 

In the continuum limit, the average particle density in the stochastic model is then given by 
(<i>A(x,t)) := ({0}|<i>^(x,t)e-^-+--^[{*->^>-{*-'^>l|{0}) 

_ / D'i'AD<S?cD^AD^c^A{^, ^)g-SA+A^c[{<I>A,c},{*A,c}] 

/ £)$^i:)$ci:)|>^£)$ce-'5-*+-*-c[{<i'A,c},{*A,c}] 

(31) 

where D denotes the measure of the fmictional integral and ^a,c is the shifted eigenvalue of the 
dual of the coherent state under the creation operator defined by <&a.c ■— ^*a c~^- Accordingly, 
all fields that incorporate shifted eigenvalues instead of the original eigenvalues will be denoted by 
a twiddle in the sequel. Note that the average (|3ip is performed taking into account the dynamics 
and the initial conditions — for a more detailed discussion see [23] . We already have incorporated 
the shifted initial state 

\^{t = 0)) = g/'i°2:(nA(0)l>A(t=0)+nc(0)*c(t=0))||Q|^^ ("32) 

in the shifted action S. The symbol Sa+a-^c represents the shifted action which is defined as 
follows 

Sa+a^c f^dt J d'^xi^A}^^ + {^c}^^ + Ha+a^c[{<^a},{^a},{'^c},{^c}], 

(33) 

with the shifted Hamiltonian Ha+a^c- The shifted action for the chemical reaction A + A — > C 
takes the form 

Sa+a^c[{'^a}, {^a}, {<i>c}, {^c}] = r dt j d^.T ( (*Af (^ - 5mA)<I>m 

- ^mQm - Am^'a/)) + '«(2$A + *i - *c)^^ 

d°x{nA{mA{t = 0) +fic(0)$c(i = 0)). (34) 

We want to untangle the quadratic term A linear expression in $^ can be obtained by means 
of a Gaussian transformation 
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where V[ri\A+A^c is the probabiUty distribution for a white noise ri[x,t). The Gaussian distri- 
bution reads 

^WA+A^c = e-^/o-'^*/<^"^'7^(x,*). (36) 

Now the shifted action S is hnear in and one can easily integrate out over ^a (x, t) and 
$c(x, i) in dSI]). One obtains 

(0[$A,$c])« J D<S>AD^cD7]0[^A:<i>c]S[J'A]S[J'c]P[v]A+A-,c 

OC j Dt] 0[^A[v{^,t),K,t],^c[vi.^,t),X,t]]P[lj]A+A^C, 

(37) 

where d[J-A,c] is a functional Dirac delta distribution. In its generalised Fourier representation 
it is defined by 

(5[J^] := constant [ d\{y)ei '^y^^y^^^^^y^^y\ (38) 



with z(j/) being a multicomponent field satisfying the constraint 

T[^{y),v]^Q. (39) 
Accordingly, the functions i'yi[?7(x, i), x, i] and $c[?7(x, t), x, t] satisfy the constraints 

^a[$a(x, t), X, t] = -^^^^ + ^aA$a(x, t) - 2S$^(x, i) - Aa^a(x, i) 



+ JA + iV2K$A(x, t)?7(x, i) = 0, (40) 



^c[$c(x, t), X, t] = ^^^-^ + i?cA$c(x, i) + K$l(x, t) - Ac$c + = 0. 

at 



(41) 



It follows from equation ([57)1 that the average particle density for the A respectively C molecules 
is now given by 

($Ac(x,t)) = / i?77$A,c[?7(x,i),x,i]e-^/o^'^*/''°-'''("'*). (42) 



The stochastic noise rj has zero mean value and a correlation given by 

(77(x,i)r;(x',t'))pM.+.^c -'^(^^(x-xXt-O- (43) 

This is obvious when considering the Gaussian distribution . Note that the above average is 
no longer taken over the initial conditions. 

The constraint equation (|40p is an inhomogeneous partial stochastic differential equation with 
additive noise for a complex fluctuating unknown field in the Ito calculus. It resembles the 
deterministic partial differential equation that describes the evolution of the mean particle density 
in the classical theory. But despite the temptation for an intuitive interpretation it is very 
important to keep in mind that in equation (|40p we are confronted with a complex fluctuating 
quantity that has as such no physical interpretation. Only if the path integral average (PIA) 
(|^^ of a solution to or (I4ip is taken over all possible realisations of the stochastic noise that 
appears in the constraint equation can one interpret the outcome of this computation as a mean 
particle density. 



3 Case A: Vanishing Source Rate 

In the remainder of this paper, let us concentrate on a single spatial site model. We will now 
compare the results in the stochastic model to the observations made in the traditional approach. 
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The classical equation for the evolution of the mean particle density in the single spatial site 
model reads 

+ 2fin\{t) + XAUAit) = (44) 

where n^(t) denotes the mean particle density of the A molecules in the mean field approximation. 
The classical evolution equation (144]) is solved by 

UA (t) = — ^-^ . (45) 

2k _i + eA^*(l + |An^(0)) 

The stochastic constraint equation for the complex fluctuating field ^A{t) in zero spatial dimen- 
sions with vanishing source rate takes the form 

^l'A(t) + 2K^l{t) + XAUAit) - iV2R^A{tHt) = 0. (46) 

The traditional equation for the average particle density of the A molecules in the mean field 
approach (|44p and the stochastic constraint equation associated with the A molecules (|46p resem- 
ble each other at first sight. But as mentioned before the solution of the stochastic differential 
equation (|46p is a complex, fluctuating field that can only be interpreted as an average particle 
density after it has been averaged in the sense of equation (|42|) . For vanishing source rate it is 
fairly easy to find an analytic solution of equation (j46|) . The stochastic constraint equation (|46|) 
— because of the continuous but not smooth nature of a stochastic process — has to be understood 
in terms of a stochastic integral equation 

^A{t) - ^A{t = 0)= f ds a($A(s), s)+ f ds b{^A{s),s)r){s), (47) 
Jo Jo 

where 

a{^A{t), t) = Ja - XA^A{t) - 2R,^\{t) : drift coefficient, 



b{^Ait),t) = W2K^A{t) : diffusion coefficient. (48) 

The stochastic noise ri{t) is rewritten in terms of the Wiener process W(t) 

ri{t)dt = dW{t). (49) 

For vanishing source rate ja the stochastic constraint equation for the A particle density ((iD|) 
reduces to the following equation 

d^A{t) = ( - 2R^\{t) - XA^A{t))dt + iV2R^ A{t)dW{t). (50) 

The above equation is a nonlinear reducible stochastic differential equation with polynomial drift 
of degree two in the Ito picture. In contrary to a Stratonovich stochastic differential equation an 
Ito stochastic differential equation can not be solved directly by methods of classical calculufl 
For an analytical solution of an Ito stochastic differential equation one has to use a modified 
version of the drift coefficient 

a{^A{t),t) — > a($A(i),t) - h(^A{t),t)^J-j-^b(^A{t),t), (51) 

where the derivative in the last term is a functional derivative. Equation (jSOp is a stochastic 
version of a Verhulst-like equation — see [M] . It can be reduced to a linear stochastic differential 
equation with multiplicative noise. We obtain the solution to the first stochastic constraint 
equation (|40|) for vanishing source rate, namely 

^A{t) = ' , . ^ . (52) 

1 + 2k$a(0) /o e(«-^-*)«+iv^w/^(^)ds 



^Sample paths of a Wiener process are — with reasonable certainty — neither differentiable nor of bounded 
variation. As a consequence one is left with different interpretations of stochastic equations, namely the Ito and 
the Stratonovich interpretation. For a further reading we refer to |14) . 
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Inserting the above solution into the path integral average (j42|l one obtains the average particle 
density for the A molecules in the stochastic picture. 

The stochastic constraint equation for the reaction product, the C particles, in zero dimen- 
sions looks — in its form — identical to the classical evolution equation for the mean density of C 
particles 

^nc (t) + Xcnc (t) - Rn\ (t) -jc = (53) 

In the single spatial site model, the full solution of the second constraint equation — simply take 
$c(i) instead of nc(t) in the above equation (|53p — can be obtained even for non- vanishing source 
rate and reads ^ 

^c{t)= (^^ e^^^(s$^(s)+jc)ds + $c(0))e-^'^*. (54) 

The stochasticity of the above equation is hidden in the first term employing the fluctuating 
solution $yi(i) of the first constraint equation (f46|) . 

Once the solutions to the stochastic constraint equations (|52p and ([5^ are known, one has to 
insert either of the solutions into the path integral average (j42p and compute the path integral 
by means of a Monte Carlo calculation in order to obtain the average particle density for the A 
or C particle population, respectively. Random samples are generated according to the Gaussian 
probability distribution ((36|) : that is we generate Wiener processes. We estimate the path integral 
(|42p by summing a large number of solutions of the constraint equations associated to the set 
of random samples generated in the above sense and divide the sum by the number of random 
samples. The Monte Carlo method displays a convergence of where N is the number of 
random samples — see |19) . 

Instead of using the expressions of the analytical solutions, equation ([5^ and equation ([5^ . 
to generate solutions to the stochastic constraint equations one can alternatively compute the 
paths directly from the stochastic differential equations (|40|) and (|4T|l . The latter method turns 
out to be less time consuming. The stochastic differential equation (j40p in zero dimensions can 
be converted into 

XA,n+l = Xa^u + (-Aa^A,™ - 2SXi „)A„ + iV2RXA,nAWn, (55) 

where Xa^u ■— ^A{tn) in discretised time tn for n — 0, ..,N, A„ := tn+i — tn and /SWn ■— W^t„+i — 
Wt^. The AW„ is generated by two uniformly distributed independently random variables via the 
Box-MuUer transformation — see, for example, [14 . The numerical scheme (l55t is called the Euler 
scheme and is the most straightforward approach to undertake some numerical investigations. 
Accordingly, the second constraint equation in the single spatial site model and for vanishing 
source rate takes the following form 

Xc,n+1 = Xc,n + i~\cXc,n + kX^,„)A„. (56) 

As stochastic differential equations are extremely sensitive one has to convince oneself that the 
code is stable and converging as it should be. Other schemes we used that might be more accurate 
or stable than the Euler method are the Milstcin scheme, the simplified order 2.0 weak Taylor 
scheme, the implicit order 1.0 strong Runge-Kutta scheme, the predictor-corrector method of 
order 1 with modified trapezoidal method weak order 1.0 — see [13]. 

For the numerical evaluation we employ values for the rate coefficients that can be found in 
realistic physical set-ups. Instead of using the coefficients introduced in the continuum limit, 
namely Xa,c smd k, we employ the traditional rate coefficients which we denote by La,c and 
Ka,c and which have dimensions per unit time. Dimensional analysis shows that using these 
new constants we are now calculating an average particle population instead of an average par- 
ticle density — this can be easily verified in equations (|46|) and (|53p. 
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The following plots were generated for the situation where two hydrogen atoms react on the 
surface of an interstellar dust particle. According to jS] [22] the reaction rate takes the value of 
K = 1.45 X 10^ s~^, the evaporation rate for the hydrogen atoms Lh — 1.88 x 10~^ and 
the evaporation rate for the reaction product — 6.9 x 10^^ s^^ . As initial values we used 
$a(0) = $c(0) = 6. 

In Figure 1 we generate the real part of one solution to the stochastic constraint equation for 
the hydrogen atoms under the above conditions. Figure 2 shows the imaginary part of the same 
solution of the stochastic constraint equation for the reaction partners. If one compares these 
plots to Figure 3 and Figure 4 where the real and imaginary part of the path integral average 
over 1000 realisations of the white Gaussian noise for the H atoms are given one observes that 
the real part of the path integral average smoothes out and the fluctuations in the imaginary part 
decrease in intensity. For increasing number of paths employed in the path integral average the 
imaginary part of the PIA tends to zero. Therefore, it is safe to interpret the real part of the 
path integral average as the average particle population. 

Figures 5 and 6 show the real and imaginary part of a solution to the second stochastic equation 
that constrains the reaction products 7?2- Again the fluctuations are smoothed out in Figure 7 
and Figure 8 when the path integral average for the diatomic hydrogen is taken over 1000 reali- 
sations of the stochastic noise. 



Together with Figure 7 one can interpret the path integral average in Figure 3 in the follow- 
ing way: according to Figure 7 the chemical reaction stops after certain transient processes. The 
average population of the diatomic hydrogen is constant. The intuitive physical reason for this 
is because all the potential reaction partners H have already been used to form the reaction 
product i?2- On the other hand, from Figure 3 one sees how the plot for the hydrogen atoms 
approaches asymptotically the value one half which could be interpreted as a state consisting of 
either zero or one particle. The value of oo) = 1/2 is the lowest possible eigenvalue of the 

coherent states under the annihilation operator once the system has reaches its equilibrium. The 
rate coefficients K and La,c have to be understood in a probabilistic sense, in a similar fashion 
as it is done with, say, the mean life expectancy of a radioactive isotope. For the specific values 
we used in our calculations the reaction rate dominates over the evaporation rate. 

The discussion above can be compared with the results obtained from the solution to the classical 
evolution equation (|45p which predicts an asymptotic value n^(t) — *■ as i — > oo for La 0. 
That is, the classical model predicts the extinction of all the reactants. 



4 Case B: non- vanishing source rate 

As in the previous section, we generate solutions to the constraint equations (I40p and (|4ip in zero 
dimensions by the numerical schemes discussed above but now with a non-zero source rate. We 
compare the results to the solutions of the classical evolution equations which read 

-AA + atanh(fft + /3)) 
nA{t) = , (57) 

where a :— \J SkJa + X\ and (3 := 2/ a arctanh((4KnA(0) + \A)/a) and 

ncW = e~^^*( / e'^^%jc + Kn\{s))ds + nc{0)). (58) 



We now analyse the influence of the source rate coefficient associated with the adsorption of 
reactants onto the surface of the seed particle. For convenience, we will leave the source rate 
for the reaction products at zero as it will not have significant influence on the outcome of our 
discussion. The other rate coefficients were chosen as before, Ja,c, La,c and Ka,c denoting the 
rate coefficients per unit time. The plots in Figures 9 to 12 were obtained for a source rate that 
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Figure 1: The real part of one possible solution to the first constraint equation (f40|) for the 
hydrogen atoms under interstellar space conditions {K = 1.45 x lO^s^^, Lh = 1-88 x 10^'^s^^) 
with vanishing source rates {Jh — Jh2 = Os""'^). 
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Figure 2: The imaginary part of the solution to the first constraint equation (|40l) for the same 
stochastic noise for the reaction partners H under interstellar space conditions {K = 1.45x 10^s~^, 
Lh = 1.88 X 10~^s~^) and for zero source rate [Jh = Jh2 = Os~^). 



12 




1e-05 2e-05 



3e-05 4e-05 5e-05 
time in seconds 



6e-05 7e-05 8e-05 



Figure 3: The real part of the path integral average (PI A) of solutions to the first constraint 
equation (|40|l with vanishing source rate {Jh — Jh2 = Os^^) for the hydrogen atoms under 



interstellar conditions (K — 1.45 x 10 s , Lh = 1- 



10- 



over 1000 possible paths. 
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Figure 4: The imaginary part of the path integral average (PIA) of solutions to the first constraint 
equation (|40p with zero source rate {Jh = Jh^ = Os^^) for the hydrogen atoms under interstellar 
conditions [K = 1.45 x IQ^s'^, Lh = 1.88 x lO-^g-i) over 1000 paths. 
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Figure 5: The real part of a possible solution to the second constraint equation (I41|l for vanishing 
source rate {Jh = Jh2 ~ Os^^), that is the constraint equation for the diatomic hydrogen under 
interstellar conditions (A' — 1.45 x 10^s^\ = 6.9 x lO^^s^^). 
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Figure 6: The imaginary part of a possible solution for one specific realisation of the stochastic 
noise to the second constraint equation (|41[) with zero source rate {Jh = Jh2 = Os^^) under 
interstellar conditions {K = 1.45 x 10^s^\ = 6.9 x 10~^s~^). 
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Figure 7: The real part of the path integral average (PIA) of solutions to the second constraint 
equation (|4T|l for the reaction product H2 under interstellar conditions {K = 1.45 x 10^s~^, 
= 6-9 X 10~^s~^) over 1000 paths with vanishing source rate (Jh — Jh2 = Os^^). 
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Figure 8: The imaginary part of the path integral average (PIA) of solutions to the second 
constraint equation ((4T|) for the diatomic hydrogen under interstellar conditions {K — 1.45 x 
lO^s^^, = 6.9 X lO^^s"^) over 1000 realisations of the stochastic noise for zero source rate 
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is big in comparison to the other rate coefhcients, Jh — 10^ s^^, whereas in Figures 13 to 16 the 
source rate was chosen to be small, Jh — 10"^ s~^. Although one observes fluctuations both in 
the real part of a single solution to the first constraint equation (Figure 9) and in the real part of 
one path associated with the second constraint equation (Figure 11), Figures 9 to 12 reproduce 
deterministic behaviour. The real part of the path integral average of the reactants (Figure 10), 
as well as the real part of the path integral average of the reaction products (Figure 12) over 1000 
possible realisations of the stochastic noise, coincide with the results of the associated classical 
equations. This accordance between classical and stochastic model is no longer valid for Figures 
13 to 16. As an example, we generated one single path for each particle population, the reaction 
partners and the reaction products, and plotted their real part in Figure 13 and Figure 15, re- 
spectively. 

Let us now compare the average particle populations of the hydrogen atoms and the diatomic 
hydrogen for large source rate (Figure 10 and Figure 12) to the average particle populations of 
the reactants and the reaction product for small source rate (Figure 14 and Figure 16). In the 
deterministic case, that is for large source rate, the chemical reaction does not die out after a cer- 
tain period of time — see Figure 12 — in contrary to the observations made from Figure 16 where 
the real part of the path integral average of the diatomic hydrogen over 1000 paths is plotted for 
a small source rate. As can be seen from Figure 10 and Figure 14 respectively, for a source rate 
of Jh = 10^ the average particle population of the H atoms reaches an asymptotic value of 
22.28 whereas for a source rate of Jh = 10~^ the average particle population of the H atoms 
is 1/2 as was the case for vanishing source rate in the latter section. Figure 12 and Figure 16 
give the average particle population for the H2 atoms at time t = 8 x 10"^ s~^, namely 3996.89 
for large source rate Jh = 10® and 8.73 for small source rate Jh = 10^^ s~^. To determine 
the transition from a deterministic to a stochastic behaviour we computed the average particle 
population of the reaction partners and products for a source rate of the hydrogen atoms in the 
range of Jh S [10® s~^, 10"'' s~^] for each order of magnitude and generated Figure 17 to Figure 
19. Already when the reaction rate and the source rate are of the same order of magnitude one 
can observe deviations from the classical behaviour, in the sense that the equilibrium value of 
the average particle population obtained by the stochastic methods is higher than that predicted 
from the classical equation ((57)) . One also observes from Figure 17 that for a source rate between 
Jh = 10^ and Jh = 10^ s""'^ the chemical reaction dies out because there is not one pair of 
reaction partners left in order to initiate a chemical reaction. 

5 Conclusions 

We have argued that the classical evolution equations for the mean particle population of a 
chemical species involved in a heterogeneous chemical reaction do not give the right results for 
small systems. Instead, we developed a stochastic model that includes statistical fluctuations and 
showed in our numerical investigations that those fluctuations can not be ignored for low rates of 
particle adsorption onto the surface of a grain particle. Although one starts from an apparently 
classical system the introduction of a quantum field theoretical formalism for its description seems 
to force us to adopt a "quantum mechanical- like" interpretation of the results. 

It is possible to extend this work to other chemical reactions, for example of the type A+B — > C. 
In a next step, we shall consider a network of chemical reactions in which several reactions com- 
pete against each other. One would expect that it will take considerably longer to reach an 
asymptotic value for the average particle population of a certain species. 
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Figure 9: Real part of one possible solution to the first constraint equation (HDl) for a value of the 
source rate of Jr = 10* for the reactants and for K = 1.45 x 10^s~^, Lh — 1-88 x 10~'^s~^. 
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Figure 10: Real part of the path integral average (PIA) of solutions to the constraint equation 
(j40l) with Jh = 10* s-^ for the hydrogen atoms and K ^ 1.45 x 10^s"\ Lh = 1.88 x IQ-^s"^ 
averaged over 1000 paths. 
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Figure 11: Real part of one path for the constraint equation (f4T|) associated with the reaction 
product i/2 plotted for the source rate Jh = 10^ s^^ and K = 1.45 x 10^s^\ Lj^^ = 6.9 x lO^^s"^, 
Jh, = Os-\ 
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Figure 12: Real part of the path integral average (PIA) for the population of diatomic hydrogen 
for a large source rate compared to the other rate coefficients, namely for Jh = 10^ s^^, and for 
K = 1.45 X 10S-\ Lff, = 6.9 x IQ-^s-i, J^, = Os-\ 
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Figure 13: Real part of a solution to the stochastic equation (PD|) constraining the reaction 
partners {K = 1.45 x 10^s~^, Lh — 1-88 x 10~^s"^) with a source rate for the hydrogen atoms 
of value Jh = 10^^ s^^. 
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Figure 14: Real part of the path integral average (PIA) of solutions to the constraint equation 
(PO)) for the hydrogen atoms over 1000 paths for a source rate of = 10^^ s^^ and for K = 
1.45 X lO^s-i, Lh = 1.88 x IQ-^s-i. 
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Figure 15: Real part of one solution to the second constraint equation (HT]) for a small source 
rate for the reactants compared to the other rate coefficients: Jh = 10^^ s^^, and with K — 
1.45 X lO^s-i, Lh^ = 6.9 X 10-^?-^, Jh^ = Os-^. 
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Figure 16: Real part of the path integral average (PIA) of the reaction products with K = 
1.45 X lO^s"^, = 6.9 X 10~*s~^, = Os^^ over 1000 realisations of the stochastic noise for 
a source rate of the hydrogen atoms of Jh ~ 10^^ s~^. 
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Figure 17: Dependence of the real part of the path integral average (PIA) of solutions to the 
first constraint equation (HDl) taken over 1000 paths after the transient processes on the source 
rate of the reactants for the average population of atomic hydrogen {K = 1.45 x 10^s~^, Lh = 
1.88 X lO-^s-i). 
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Figure 18: Real part of the path integral average (PIA) of solutions to the second constraint 
equation gl|) = 1.45 x 10S-\ L//, = 6.9 x 10-*s-\ Jh, = Os'^) over 1000 paths for the 
average population of molecular hydrogen at t — 8 x 10^^ s versus the logarithm of the source 
rate of the hydrogen atoms. 
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Figure 19: Same curve as shown in Figure 18 but zoomed to resolve the plot close to the horizontal 
axis. 
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